% clear all;
% clc;

patnum = 1;
dir1=['patient',num2str(patnum),'_8_12_13'];
dir2='other_geom';
dir3=['p',num2str(patnum),'_act_curves'];
addpath(dir1);
addpath([dir1,'\',dir2]);
addpath([dir1,'\',dir3]);

% ----------------------- Importing Data ----------------------------------

noaxons=200;
epidur=0;
th_elec=-20; % 0, -10, and -20
patgeom=load('patient_cordgeom.txt');

if(epidur==1)
    loctag='epi_';
    ftag=['d',strrep(num2str(0.3),'.','p'),'_'];
else
    loctag='sub_';
    ftag=['d',strrep(num2str(0.2),'.','p'),'_'];
end
thtag=['t',strrep(num2str(th_elec),'-','n')];
dcpre=[loctag,'LT6_dc_act_p',num2str(patnum),'_',ftag,thtag];
drpre=[loctag,'LT6_dr_act_p',num2str(patnum),'_',ftag,thtag];

uni_cur=load([loctag,'atp_current_p',num2str(patnum),'_',ftag,'t0','.txt']);

dc_th=abs( load(['pos_',dcpre,'.txt']) );
dc_th(:,1)=dc_th(:,1)/(2/uni_cur)*1e3;
dcbyderm = zeros(10,10);
dcbyderm(:) = dc_th(1:(noaxons/2),1);
dr_th=abs( load(['pos_',drpre,'.txt']) );
dr_th(:,1)=dr_th(:,1)/(2/uni_cur)*1e3;

dcth_sel=sort(dc_th(:,1));
drth_sel=sort(dr_th(:,1));
dr0=100*sum(dcth_sel<drth_sel(1))/noaxons;

aa = 1:10;
bb = (aa'*ones(1,10))';
xx = bb(:);
yy = dcbyderm(:);

numax=10;
for ii=1:numax
    
    ko=numax*(ii-1)+1;
    kf=numax*ii;
    
    for jj=1:(numax-1)
        
        kk=numax*(ii-1)+jj;
        yycomp=[zeros(jj,1);yy((kk+1):kf)];

        ck1=aa(yy(kk)==yycomp);
        
        if(any(ck1))
            kk_p=numax*(ii-1)+ck1;
            xx(kk_p)=xx(kk_p)+(1:length(ck1))'*(2/numax);
        end
        
    end
end

% rtt1 = 4/0.011493/1e3;
% ratp = 4/0.005924/1e3;
% yy = yy/rtt1;

% ---------------------- Plotting Results ---------------------------------

thtag = ['Lateral Deviation = ',num2str(th_elec),'^o'];
t8_vertlv  = {'T9','T10/11','T12','L1','L2/3','L3/4','L5','S1/2','S2/3','S4/5'};

figure;
hold on;
plot(xx,yy,'k.','MarkerSize',30);
plot([0,11],[min(dr_th(:,1)),min(dr_th(:,1))],'k--','LineWidth',4);
ttl = ['Patient ',num2str(patnum),', ',thtag];
title(ttl,'FontSize',30);
ylabel('Stimulation Threshold Current (mA)','FontSize',36);
hold off;
xlim([0,11]);
ylim([0.05,1]);
set(gca,'XTick',[],'YScale','Log','YTick',[0,0.1,0.2,0.5,1,1.5,2],'FontSize',36);
%set(gca,'YScale','Log');
for jj = 1:length(t8_vertlv)
    nolet = length(t8_vertlv{jj});
    hjj=text(jj-0.175*nolet/2,0.035,t8_vertlv{jj},'FontSize',30);
    set(hjj,'rotation',60);
end
legend('DC fibers','DC_0','Location','NE');
axis square;
